Consider the following data about sets of brothers, each set consisting of an older and younger brother.
brothers <- data.frame(
older = c(21, 10, 50, 15, 46, 39, 35, 25, 39, 31, 40, 35, 45, 33),
younger = c(44, 9, 61, 19, 44, 43, 26, 27, 26, 40, 57, 56, 52, 39)
)
The question we’re considering is whether or not there is a difference between older and younger brothers. One way we can quantify this difference is by taking the difference of the median value for each group.
(observed_difference <- median(brothers$older) - median(brothers$younger))
[1] -6.5
Now, how do we determine if this is “significant”? Using a permutation test, we can create a distribution of data under the null hypothesis. That is, we can simulate scenarios where the null hypothesis is true.
In this case, if the null hypothesis is true, that means that there is no difference between older and younger brothers. If that’s the case, it shouldn’t matter which brother is assigned which score. This is the basis of the permutation test:
The distribution of the test statistic under the null hypothesis is obtained by calculating all possible values of the test statistic under rearrangements of the labels on the observed data points (Wikipedia).
Essentially, we’re going to create the null hypothesis by breaking the connection between label and response.
To build some intuition behind this, let’s walk through a single permutation of the brothers data.
(perm_sample <- apply(brothers, 1, sample))
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]
[1,] 44 9 61 19 46 43 35 25 26 31 40 35 52 39
[2,] 21 10 50 15 44 39 26 27 39 40 57 56 45 33
Now that we’ve gone through a single permutation, let’s calculate the difference in medians given this version of the data.
median(perm_sample[1,]) - median(perm_sample[2,])
[1] -2
Now, following the theme of this class, we’re going to do that “a bunch of times” so that we can build up a distribution of values under the null hypothesis.
n_permutations <- 10000
results <- replicate(n_permutations, {
perm_sample <- apply(brothers, 1, sample)
median(perm_sample[1,]) - median(perm_sample[2,])
})
Now we have results which contains 10,000 calculated differences under the null hypothesis. We can look at a distribution of these values.
plot(density(results))
We can also plot our original observed_difference.
plot(density(results))
abline(v = observed_difference, col = "red")
Now, given this, how can we answer our original question (is there a difference)? Well, knowing the definition of a p value, we can do the following:
(p_value <- mean(results <= observed_difference))
[1] 0.0648
Now, are we done? Not quite! Remember, this p_value is only an estimate! We need to assess and report on the uncertainty around this estimate.
ci <- p_value + c(-1, 1) * qnorm(.975) * sqrt(p_value * (1 - p_value) / n_permutations)
c(lower = ci[1],
p_value = p_value,
upper = ci[2])
lower p_value upper
0.05997511 0.06480000 0.06962489
LS0tCnRpdGxlOiAiQnJvdGhlcnMiCm91dHB1dDogaHRtbF9ub3RlYm9vawotLS0KCmBgYHtyIHNldHVwLCBpbmNsdWRlID0gRkFMU0V9CmtuaXRyOjpvcHRzX2NodW5rJHNldChjYWNoZSA9IFRSVUUpCmBgYAoKQ29uc2lkZXIgdGhlIGZvbGxvd2luZyBkYXRhIGFib3V0IHNldHMgb2YgYnJvdGhlcnMsIGVhY2ggc2V0IGNvbnNpc3Rpbmcgb2YgYW4Kb2xkZXIgYW5kIHlvdW5nZXIgYnJvdGhlci4KCmBgYHtyfQpicm90aGVycyA8LSBkYXRhLmZyYW1lKAogIG9sZGVyID0gYygyMSwgMTAsIDUwLCAxNSwgNDYsIDM5LCAzNSwgMjUsIDM5LCAzMSwgNDAsIDM1LCA0NSwgMzMpLAogIHlvdW5nZXIgPSBjKDQ0LCA5LCA2MSwgMTksIDQ0LCA0MywgMjYsIDI3LCAyNiwgNDAsIDU3LCA1NiwgNTIsIDM5KQopCmBgYAoKVGhlIHF1ZXN0aW9uIHdlJ3JlIGNvbnNpZGVyaW5nIGlzIHdoZXRoZXIgb3Igbm90IHRoZXJlIGlzIGEgZGlmZmVyZW5jZSBiZXR3ZWVuCm9sZGVyIGFuZCB5b3VuZ2VyIGJyb3RoZXJzLiBPbmUgd2F5IHdlIGNhbiBxdWFudGlmeSB0aGlzIGRpZmZlcmVuY2UgaXMgYnkgdGFraW5nCnRoZSBkaWZmZXJlbmNlIG9mIHRoZSBtZWRpYW4gdmFsdWUgZm9yIGVhY2ggZ3JvdXAuCgpgYGB7cn0KKG9ic2VydmVkX2RpZmZlcmVuY2UgPC0gbWVkaWFuKGJyb3RoZXJzJG9sZGVyKSAtIG1lZGlhbihicm90aGVycyR5b3VuZ2VyKSkKYGBgCgpOb3csIGhvdyBkbyB3ZSBkZXRlcm1pbmUgaWYgdGhpcyBpcyAic2lnbmlmaWNhbnQiPyBVc2luZyBhIHBlcm11dGF0aW9uIHRlc3QsIHdlCmNhbiBjcmVhdGUgYSBkaXN0cmlidXRpb24gb2YgZGF0YSB1bmRlciB0aGUgbnVsbCBoeXBvdGhlc2lzLiBUaGF0IGlzLCB3ZSBjYW4Kc2ltdWxhdGUgc2NlbmFyaW9zIHdoZXJlIHRoZSBudWxsIGh5cG90aGVzaXMgaXMgdHJ1ZS4KCkluIHRoaXMgY2FzZSwgaWYgdGhlIG51bGwgaHlwb3RoZXNpcyBpcyB0cnVlLCB0aGF0IG1lYW5zIHRoYXQgdGhlcmUgaXMgbm8KZGlmZmVyZW5jZSBiZXR3ZWVuIG9sZGVyIGFuZCB5b3VuZ2VyIGJyb3RoZXJzLiBJZiB0aGF0J3MgdGhlIGNhc2UsIGl0IHNob3VsZG4ndAptYXR0ZXIgd2hpY2ggYnJvdGhlciBpcyBhc3NpZ25lZCB3aGljaCBzY29yZS4gKipUaGlzIGlzIHRoZSBiYXNpcyBvZiB0aGUKcGVybXV0YXRpb24gdGVzdDoqKgoKPiBUaGUgZGlzdHJpYnV0aW9uIG9mIHRoZSB0ZXN0IHN0YXRpc3RpYyB1bmRlciB0aGUgbnVsbCBoeXBvdGhlc2lzIGlzIG9idGFpbmVkIApieSBjYWxjdWxhdGluZyBhbGwgcG9zc2libGUgdmFsdWVzIG9mIHRoZSB0ZXN0IHN0YXRpc3RpYyB1bmRlciByZWFycmFuZ2VtZW50cyBvZiAKdGhlIGxhYmVscyBvbiB0aGUgb2JzZXJ2ZWQgZGF0YSBwb2ludHMgKFtXaWtpcGVkaWFdKGh0dHBzOi8vZW4ud2lraXBlZGlhLm9yZy93aWtpL1Jlc2FtcGxpbmdfKHN0YXRpc3RpY3MpI1Blcm11dGF0aW9uX3Rlc3RzKSkuCgoKRXNzZW50aWFsbHksIHdlJ3JlIGdvaW5nIHRvIGNyZWF0ZSB0aGUgbnVsbCBoeXBvdGhlc2lzIGJ5IGJyZWFraW5nIHRoZQpjb25uZWN0aW9uIGJldHdlZW4gbGFiZWwgYW5kIHJlc3BvbnNlLgoKVG8gYnVpbGQgc29tZSBpbnR1aXRpb24gYmVoaW5kIHRoaXMsIGxldCdzIHdhbGsgdGhyb3VnaCBhIHNpbmdsZSBwZXJtdXRhdGlvbiBvZiAKdGhlIGBicm90aGVyc2AgZGF0YS4KCmBgYHtyfQoocGVybV9zYW1wbGUgPC0gYXBwbHkoYnJvdGhlcnMsIDEsIHNhbXBsZSkpCmBgYAoKTm93IHRoYXQgd2UndmUgZ29uZSB0aHJvdWdoIGEgc2luZ2xlIHBlcm11dGF0aW9uLCBsZXQncyBjYWxjdWxhdGUgdGhlIGRpZmZlcmVuY2UKaW4gbWVkaWFucyBnaXZlbiB0aGlzIHZlcnNpb24gb2YgdGhlIGRhdGEuCgpgYGB7cn0KbWVkaWFuKHBlcm1fc2FtcGxlWzEsXSkgLSBtZWRpYW4ocGVybV9zYW1wbGVbMixdKQpgYGAKCk5vdywgZm9sbG93aW5nIHRoZSB0aGVtZSBvZiB0aGlzIGNsYXNzLCB3ZSdyZSBnb2luZyB0byBkbyB0aGF0ICJhIGJ1bmNoIG9mCnRpbWVzIiBzbyB0aGF0IHdlIGNhbiBidWlsZCB1cCBhIGRpc3RyaWJ1dGlvbiBvZiB2YWx1ZXMgdW5kZXIgdGhlIG51bGwKaHlwb3RoZXNpcy4KCmBgYHtyfQpuX3Blcm11dGF0aW9ucyA8LSAxMDAwMApyZXN1bHRzIDwtIHJlcGxpY2F0ZShuX3Blcm11dGF0aW9ucywgewogIHBlcm1fc2FtcGxlIDwtIGFwcGx5KGJyb3RoZXJzLCAxLCBzYW1wbGUpCiAgbWVkaWFuKHBlcm1fc2FtcGxlWzEsXSkgLSBtZWRpYW4ocGVybV9zYW1wbGVbMixdKQp9KQpgYGAKCk5vdyB3ZSBoYXZlIGByZXN1bHRzYCB3aGljaCBjb250YWlucyAxMCwwMDAgY2FsY3VsYXRlZCBkaWZmZXJlbmNlcyAqdW5kZXIgdGhlCm51bGwgaHlwb3RoZXNpcyouIFdlIGNhbiBsb29rIGF0IGEgZGlzdHJpYnV0aW9uIG9mIHRoZXNlIHZhbHVlcy4KCmBgYHtyfQpwbG90KGRlbnNpdHkocmVzdWx0cykpCmBgYAoKV2UgY2FuIGFsc28gcGxvdCBvdXIgb3JpZ2luYWwgYG9ic2VydmVkX2RpZmZlcmVuY2VgLgoKYGBge3J9CnBsb3QoZGVuc2l0eShyZXN1bHRzKSkKYWJsaW5lKHYgPSBvYnNlcnZlZF9kaWZmZXJlbmNlLCBjb2wgPSAicmVkIikKYGBgCgpOb3csIGdpdmVuIHRoaXMsIGhvdyBjYW4gd2UgYW5zd2VyIG91ciBvcmlnaW5hbCBxdWVzdGlvbiAoaXMgdGhlcmUgYQpkaWZmZXJlbmNlKT8KV2VsbCwga25vd2luZyB0aGUgZGVmaW5pdGlvbiBvZiBhIHAgdmFsdWUsIHdlIGNhbiBkbyB0aGUgZm9sbG93aW5nOgoKYGBge3J9CihwX3ZhbHVlIDwtIG1lYW4ocmVzdWx0cyA8PSBvYnNlcnZlZF9kaWZmZXJlbmNlKSkKYGBgCgpOb3csIGFyZSB3ZSBkb25lPyBOb3QgcXVpdGUhIFJlbWVtYmVyLCB0aGlzIHBfdmFsdWUgaXMgb25seSBhbiBlc3RpbWF0ZSEgV2UgbmVlZAp0byBhc3Nlc3MgYW5kIHJlcG9ydCBvbiB0aGUgdW5jZXJ0YWludHkgYXJvdW5kIHRoaXMgZXN0aW1hdGUuCgpgYGB7cn0KY2kgPC0gcF92YWx1ZSArIGMoLTEsIDEpICogcW5vcm0oLjk3NSkgKiBzcXJ0KHBfdmFsdWUgKiAoMSAtIHBfdmFsdWUpIC8gbl9wZXJtdXRhdGlvbnMpCmMobG93ZXIgPSBjaVsxXSwKICBwX3ZhbHVlID0gcF92YWx1ZSwKICB1cHBlciA9IGNpWzJdKQpgYGAKCg==